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ABSTRACT 

Population synthesis studies constitute a powerful method to reconstruct the birth 
distribution of periods and magnetic fields of the pulsar population. When this method 
is applied to populations in different wavelengths, it can break the degeneracy in 
the inferred properties of initial distributions that arises from single-band studies. In 
this context, we extend previous works to include X-ray thermal emitting pulsars 
within the same evolutionary model as radio-pulsars. We find that the cumulative 
distribution of the number of X-ray pulsars can be well reproduced by several models 
that, simultaneously, reproduce the characteristics of the radio-pulsar distribution. 
However, even considering the most favourable magneto-thermal evolution models 
with fast field decay, log-normal distributions of the initial magnetic field over-predict 
the number of visible sources with periods longer than 12 s. We then show that the 
problem can be solved with different distributions of magnetic field, such as a truncated 
log-normal distribution, or a binormal distribution with two distinct populations. We 
use the observational lack of isolated NSs with spin periods P > 12 s to establish 
an upper limit to the fraction of magnetars born with B > 10 15 G (less than 1%). 
As future detections keep increasing the magnetar and high-B pulsar statistics, our 
approach can be used to establish a severe constraint on the maximum magnetic field 
at birth of NSs. 
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1 INTRODUCTION 

The last couple of decades have seen tremendous advances in 
our understanding of the properties and evolutionary paths 
of neutron stars (NSs), largely driven by multi-wavelength 
detection. Despite this progress, however, several puzzles 
still remain. In particular, the apparent diversity of the 
observational manifestations of NSs, which has led over 
the years to several subclassifications (i.e. rotation-powered 
pulsars, dim isolated NSs, Anomalous X-ray Pulsars, Soft- 
Gamma Ray Repeaters, Central Compact Objects, among 
others) has highlighted the need of understanding possible 
evolutionary links among different classes, as well as the role 
played by the birth properties of the NSs. 

Studying the properties of NSs as a whole is best done 
via Monte Carlo simulations aimed at reproducing the galac¬ 
tic population. Since the fraction of pulsars that can be de¬ 
tected close to their birth constitutes a negligible fraction of 


the total sample, population synthesis studies (PS in the fol¬ 
lowing) generally use the present-day observed properties of 
pulsars, together with some assumptions about their time 
evolution, to reconstruct the birth distribution of periods 
and magnetic fields for the pulsar population. 


PS modeling of NSs has a long history, 
and especially so in the radio band, where the 
largest sample of pul sar s has been detecte d (e.g . 
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years into this area of research stem from the fact that 
these studies allow to constrain the birth properties of 
NSs, which are in turn intimately related to the physical 
processes occurring during the supernova (SN) explosion 
and in the proto-NS. As such, they bear crucial information 
on the physics of core-collapse SNe, in which most NSs are 
believed to be formed. 

The energy reservoir powering Rotation Powered Pul¬ 
sars (RPPs) is provided by the loss of rotational energy. 
RPPs are mostly seen in the radio band (~ 2200 ob¬ 
jects), and we can detect more than one hundred of them 
also in the X-ray and/or 7 -ray frequency range. Part of 
the rotational energy loss is converted into non-thermal 
synchrotron and curvature radiation by acceleration of 
charged particles, which are either extracted from the sur¬ 
face or created by pair production. According to the classi¬ 
cal scenario, such acceleration takes place in some regions, 
called gaps, where a cascade of pairs is supported by the 
same high-energy photons emitted by the particles, and is 
thought to produce the coherent radio emission we observe. 
While radio em issio n is likely prod uced above the polar 
cap (ISturrockl Il97ll : iRuderman fc Sutherland! Il975l l , high- 
energy phot ons are thought to originate in the outer mag¬ 
netosphere (ICh eng. Ho fc Rudermanl ll 9S(J : iR.omanil 19961 
Muslimov fc Hardind 120031 ") or in the wind region (IPetril 
2 012|j. Of particu lar interest is the recen t FIDO model 
( Kalapotharakos. Harding fc Kazanas 12014), according to 
which the gamma-ray emission is produced in regions near 
the equatorial current sheet. These models can fit the ob¬ 
served values of the radio-lags and reproduce the observed 
light curve phenomenology and the GeV photon cut-off en¬ 
ergies. 

The majority of PS studies in the literature have as¬ 
sumed that the magnetic field is fixed at its value at birth, or 
have adopted simple, parametrized expressions for its decay. 
Recently, a large theoretical effort has been devoted to model 
the magnetic field evolution with magnetothermal simula¬ 
tions which follow the coupled evolution of th e tempera¬ 
ture and the magnetic field in the NS interior |Vigano_et_alj 
l2013ll . We will use results from these simulations to model 
the magnetic f ield evolution. 

lOullon et al.l (feoiA (paper I, hereafter) have revisited 
PS modeling of isolated radio-pulsars incorporating the 13- 
field evolution from the state-of-the-art magneto-thermal 
evolution models, together with recent results on the evolu¬ 
tion of the angle between the magnetic and rotational axes 
from magnetospheric simulations. In an era where the avail¬ 
ability of high energy satellites has enlarged considerably our 
understanding of the NS population, unvealing a variety of 
observational manifestations of NSs, it is then timely to in¬ 
clude these further constraints in the global PS study of the 
NS population. A first attempt in this direction, including 
part of the X-ray emitting NS population and m agnet ic field 
decay models, was performed by IPqpov et al.l (l2010li . with 
the specific goal of uncovering links among the apparently 
different NS manifestations. 

In line with recent theoretical efforts aimed at un¬ 
derstanding the various observation al manifestations of 
NSs, with a u nifying evolution model ([Perna fc PondfioTi] : 

I Pons fc Pernal 1201 ll : IVigano et al.l 12013) ). our goal in this 
work is to combine previous studies of the radio-pulsar pop¬ 
ulation (paper I) with the study of the population of X-ray 


pulsars. For young and middle-age NSs (< 1 Myr), the ther¬ 
mal emission comes from the internal heat which is gradu¬ 
ally released or, in the case of NSs with large magnetic fields 
(B > 10 14 G), by Joule dissipation of the electrical currents 
circulating in the NS crust which enhances the thermal lumi¬ 
nosity. We do not include in this study non-thermal emission 
powered by losses of rotational energy, mainly because of the 
lack of detailed models for the physical mechanisms. Pre¬ 
dictions for the non-thermal luminosities are based on phe¬ 
nomenological correlations with the rotational energy losses, 
which include more free parameters that are not linked to 
physical properties. On the contrary, independent, physi¬ 
cally motivated cooling models provide the information link¬ 
ing the model parameters with the thermal emission of NSs. 
We hence focus on this population of NSs and leave the study 
of rotation powered, non-thermal emission (X-ray and 7 -ray 
pulsars) for future works. 

The article is organized as follows: In §[2]we summarise 
the main results of paper I about the radio-pulsar popula¬ 
tion and we describe the assumptions made to model the 
X-ray thermal radiation. In §[3]we discuss the results of our 
analysis and finally, in § [5] we present our conclusions and 
final remarks. 


2 MULTI-WAVELENGTH POPULATION 

SYNTHESIS 

From an observational point of view, there is an impor¬ 
tant qualitative difference, regarding detectability biases, 
between the radio band and the high-energy bands. In the 
former, relatively complete surveys allow to quantitatively 
implement the instrumental detectability in the PS study. 
In X-rays, instead, the sample is not complete, and the sky 
coverage is inhomogeneous. Moreover, the detectability de¬ 
pends on non-trivial issues, such as the presence of a coun¬ 
terpart with known ephemerides in other wavelengths, or 
whether or not it has been associated to a transient X-ray 
outburst. This implies that the potential number of visible 
NSs above a certain threshold flux should be larger than the 
actual number of detected sources. 

In X'-rays we have to account for the effect of absorption 
by the interstellar medium. The hydrogen column density, 
Nh, for each synthetic pulsar is: 

Nh = j riH dl (1) 

where uh is the local hydrogen number density, wh o se spa - 
tial distribution is estimated following IZane et akl )l995h . 
and the integral is performed along the line of sight l. 
The energ y-dependent absorpt io n is calc ulated following the 
models of iBalucinska-Church fc McCammonl (Il992h . which 
include photoelectric cross sections of 17 elements in the 
[0.1 —10] keV range. This allows us to calculate the detected 
(absorbed) flux on Earth Sx' s , assuming that radiation is 
isotropically distributed. 


2.1 The radio pulsar population 

For completeness, we begin our discussion by reviewing the 
main results of paper I, to which we refer for technical de- 
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Model 

Decay 

Envelope 

L 1 Bo 

log B [G] 

a B 0 

log B [G] 

rp 0 

M 

°Po 

[s] 

a 

?7-br 

[century -1 ] 

V 

A 

No decay 

Heavy elements 

12.65 

0.50 

0.38 

0.35 

0.50 

5.70 ±0.16 

0.087 ±0.017 

B 

Slow 

Heavy elements 

13.04 

0.55 

0.23 

0.32 

0.44 

2.63 ±0.05 

0.082 ±0.010 

C 

Fast 

Heavy elements 

13.20 

0.72 

0.37 

0.33 

0.41 

3.82 ±0.14 

0.106 ±0.011 

D 

Slow 

Light elements 

12.99 

0.56 

0.16 

0.31 

0.43 

2.63 ± 0.06 

0.078 ± 0.009 

E 

Medium 

Light elements 

13.13 

0.68 

0.32 

0.19 

0.44 

3.16 ±0.11 

0.083 ± 0.007 

F 

Slow (with toroidal field) 

Light elements 

12.99 

0.59 

0.21 

0.32 

0.44 

2.95 ± 0.09 

0.083 ± 0.008 


Table 1. Optimal parameters for the radio-pulsar population for each magnetothermal evolution model considered in this work. Columns 
4 to 8 indicate the values of the free parameters that minimize the T >-value computed through a 2D KS test. Column 9 shows the 
corresponding birthrate, and column 10 shows the average P-value of 10 realizations and its standard deviation <r. 
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Figure 1. Contour plot of the P-value computed through a 2D 
KS-test in the /ig 0 - ag 0 plane (a 2D cut of the general 5D 
parameter space), for models B (top), C (middle) and D (bottom). 
The purple crosses indicate the minima found by the simulated 
annealing procedure in each case, and listed in Table [T] 


tails. An outline of the procedure to generate a NS popula¬ 
tion is sketched in the following way: 

• A large number of NSs with uniformly distributed 
random ages are generated in the framework of a given 
magneto-thermal evolution model (we explore six different 
models), which determines the evolution of the magnetic 
field. The birth location distribution, kick velocity, and evo¬ 
lution in the galactic gravitational potential are described 
in paper I. 

• Gaussian distributions are assumed for both the initial 
spin period Po and the logarithm of the initial magnetic 
field, Bo. The central values (ub 0 , UPq) an( i widths ( ag 0 , 
(jp 0 ) are left as free parameters. The radio-luminosity (at 
1400 MHz) obeys the following equation: 

Lmoo = L 0 10 w (P" 3 P)“ , (2) 

where P is the spin period, P is the period time derivative, 
Lo = 5.69 x 10 6 mjy kpc 2 and L corr is a random correc¬ 
tion c hosen from a zero-centere d gaussian of a = 0.8 (as 
in lFaucher-Giguere fe Kaspill2006l) . to take into account the 
large observed dispersion. The index a is a free parameter. 

• The NSs are classified as radio-visible according to the 
same flux and beam visibility criteria of the radio survey 
used in paper I. The generation of pulsars stops when the 
number of radio-visible pulsars is equal to the number of 
radio pulsars in the considered observational sample. This 
fixes our normalisation for the total number of sources or, 
equivalently, the birth rate rib r for each synthetic population. 

• The set of free parameters (ub 0 , o-g 0 , fip 0 , ap 0 and 
a) is determined for each model by minimizing the associ¬ 
ated P-value of a generalized two-dimensional Kolmogorov- 
Smirnov test. The V statistic is calculated from the four nat¬ 
ural quadrants that define each point in the P — P plane of 
the observed sample. The fraction of data from both samples 
(observed and simulated) in the radio band, is calculated on 
each quadrant, and the P-value is defined as the maximum 
difference between observed and simulated fractions over all 
points and quadra nts ( for a more tec hnical de script ion we 
refer the reader to iFasano fe Franceschini|[l987l : IPress et al.l 
Il993h . 

The results for the 6 evolutionary models studied are 
summarized in Table |T| Here we considered three models 
with iron envelopes (A, B and C, already presented in paper 
I), differing in the timescale of the magnetic field decay. We 
have also used three new models (D, E and F) with light 
element envelopes. The last one also has an additional strong 
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toroidal magnetic field. Although the envelope composition 
barely affects the evolution of rotational properties (P and 
P), it has a relevant impact on the thermal luminosity (see 
next section). 

In general, we find that the differences in the initial 
period distribution between different models are not signifi¬ 
cant. Indeed, any broad distribution of initial periods in the 
range 0 < Po < 0.5 s can be reconciled with the data, but 
narrow distributions peaked at short spin periods (Po < 100 
ms) are ruled out. The V —value is more dependent on the 
initial magnetic field distribution, but this is not strongly 
constrained because its central value and width are corre¬ 
lated showing an interval of possible solutions: the larger the 
central value, the wider the distribution. This holds true for 
models with and without field decay, and with different as¬ 
sumptions for the magnetospheric torque. This correlation 
is visible in Fig. |T] where we show contour plots of the 11- 
value in the p,B 0 — cjb 0 plane. We note this is a 2D plot of 
a section of a 5D parameter space. For conciseness, we only 
show plots for models B,C, and D (other models look very 
similar). 

The dark blue regions, where the D-value is low, de¬ 
fine the range of allowed parameters. The purple crosses 
mark the different solutions for which explicit values for the 
parameters have been given in Table [T] The main conclu¬ 
sion from the analysis of only the radio pulsar population is 
that, independently of the underlying physical model, one 
can find acceptable parameter sets to fit the radio popula¬ 
tion reasonably well. Thus the radio-pulsar population alone 
cannot help much to discriminate among different evolution¬ 
ary models. The next question is whether or not the analysis 
of pulsar distributions in other bands can break the degen¬ 
eracy. 


2.2 Thermally emitting NSs 

2.2.1 The observational sample 

We consider all NSs for which the thermal emission is domi¬ 
nant, or clearly contributing to the total power in the soft X- 
ray band, and whose origin is attributed to residual cooling 
fr om bi rth. The sele cted sample we use has been described 
in lVigano et al.l (120131 1 and it is periodically updated in our 
websitcQ. includes thirteen closeby RPPs (only those with 
a clear thermal component), seven X-ray dim Isolated NSs 
(XINSs), and seventeen magnetars. 

The XINSs are a small group of nearby (within about 
200-300 parsecs), thermally emitting, isolated NSs. They are 
radio quiet but relatively bright in the X-ray band, and five 
of them have periods similar to magnetars. Moreover, their 
inferred magnetic fields are between the upper end of the 
radio pulsar population and the magnetar candidates. At 
variance with the majority of X-ray emitting pulsars, XINSs 
spectra are rather perfect blackbodies (kT~0.1 keV, occa¬ 
sionally with some broad lines), which make them the key 
objects for research on the NS equation of state, and the 
cooling properties of th e NS crust and core. 

The “magnetars” (lOlausen fe KaspHl2014l) are a small 
group of X-ray pulsars with spin periods between 0.3 — 12s, 


www.neutronstarcooling.info 


whose properties cannot be explained within the com¬ 
mon scenario for pulsars. Their inferred magnetic fields 
are B ~ 10 14 ~ 15 G. The decay and instabilities of such 
l arge ma gnetic fields are tho ught to power their emissi on 
dDuncan fe Thompsonl Il992l : IThompson fe Duncanl Il993l ) . 
Their powerful persistent X-ray emission (much larger 
than the rotational energy loss) is well modelled by a 
blackbody with kT ~ 0.2 — 0.8 keV with an additional 
power-law component (F ~ 2 — 4), usually interpreted 
as the product of Resonant Compton Scattering (RCS) 
of th e seed t hermal photons in the de nse magnetosphere 
(IThompson. Lvutikov fc Kulkariiil l2002lf . Hence even this 
apparent non-thermal component has a thermal origin. Dif¬ 
ferent physical models have succeeded to reproduce the 
effect of the RCS to fit the magnetars X-ray spectra 
dRea et a! 1 l2008l : IZane et al.l2009l : lBeloborodovll2013| , l. Given 
the transient nature of these objects (they show flares and 
sometimes long radiative X-ray outbursts), we have used in 
this work only X-ray luminosities derived during their qui¬ 
escent state (available for the 17 selected objects). 


2.2.2 Magneto-thermal evolution and emission model 


In NSs endowed with strong magnetic fields, the evolution 
of the internal and surface temperatures and the magnetic 
field are closely finked. The thermal luminosity is strongly 
dependent on the time-dependent magnetic field strength 
and topology, while the evolution of the magnetic field de¬ 
pends on the local values of the magnetic diffusivity and 
the electron rela xati o n tim e , whic h strongly depend on the 
temperature. IViaano et al.l (l2013l l performed 2D magneto- 
thermal simulations of NSs, and studied the evolution of 
magnetic field strength and thermal luminosity as a func¬ 
tion of time. Using this code, we produce a number of ta¬ 
bles with the effective temperature and value of the dipolar 
magnetic field as a function of the age of the NS. Among 
the many model parameters (initial magnetic field topology, 
mass, radius, envelope composition, etc), we will focus on 
the parameters that more strongly influence the observable 
quantities (luminosity, P, P) at a given age. These are the 
envelope composition (light vs. heavy elements), the initial 
magnetic field strength at the pole, Bo, and the impurity 
parameter in the inner-crust Ql^p, a measure of the aver¬ 
age cha rge distrib ution o_fthe nuclei present in the crust. We 
refer to lVigano et al.l |2013l ) for more details. 

The magnetic field strength determines the rotational 
evolution of the NS. In addition, in highly magnetised NSs, 
the additional energy reservoir provided by the magnetic 
field dissipation may further increase the surface tempera¬ 
ture, and extend the duration of the stage in which neutron 
stars are bright enough to be seen as X-ray pulsars. Q\m P was 
found to be a crucial parameter to understand the cluster¬ 
ing of spin periods in isolated X-ray pulsars and why there 
is an upper limit of isolat ed X-ray pulsars at about 12 s 
(IPons. Vigano fc Real l2013l 'l. The composition of the enve¬ 
lope is also of particular relevance. Generally speaking, low- 
fi NSs with light element envelopes have significantly larger 
temperatures and luminosities during the first ~ 10 s yr 
of their life, and lower temperature/lumi nosity afterwards , 
once they enter the phot o n cooling era jP agc ct al. 200J; 
Yatovley_j^j)ethiclj_ [2004l : IPotekhin, Chabrier fe Yakovlev! 
20071: IPage et al.ll201lli. 
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Figure 2. Left: Blackbody spectrum with kT = 1 keV (dashed line), and reprocessed spectra using the RCS model with fi r r = 0.3 and 
To = 0.1,0.5,1,3,5,10 (solid color lines). Unabsorbed fluxes with a typical galactic value of Njj ~ 3 X 10 22 cm -2 are plotted. Middle: 
Ratio between the absorbed flux of pure blackbody model to that of the RCS model as a function of the blackbody flux Sbb ■ Right: 
same as the middle panel, as a function of the logarithm of the magnetic field. The results show only the NSs with detectable thermal 
emission in the 0.1-10 keV band for model E, with fluxes above Sxmin = 10 -14 erg s —1 cm -2 . 


These three parameters characterize the evolution of the 
magnetic field strength and temperature as the star ages, 
and therefore the trajectories that NSs follow in the P-P di¬ 
agram. Having all these tabulated quantities computed via 
simulations, for each star generated in the synthetic sample 
we calculate the present value of magnetic field, P , P, tem¬ 
perature and luminosity. We assume blackbody emission, 
and integrate over the whole star surface to calculate the 
luminosity. For NSs with high magnetic fields, we also add 
the effect of resonant Compton scattering (see next subsec¬ 
tion) in the spectrum. Among all NSs generated (not only 
those seen as radio-pulsars), we select those for which the 
absorbed X-ray flux at the Earth location is above 1CP 15 erg 
s' 1 cm -2 , as potentially visible X-ray pulsars (for present 
or future instruments). 

2.2.3 Resonant Compton Scattering 

The thermal emission is usually assumed to be blackbody¬ 
like. However, it is well known that magnetar spectra show 
a non-thermal tail, whose origin cannot be attributed to 
the rotational energy of the star. The most popular inter¬ 
pretation is that magnetospheric plasma boosts up thermal 
photons coming from the surface. For these reasons, it is 
important to account for the possible effects that spectral 
distortions in the form of a hard tail would have on the 
detectability of NSs in X-rays. 

We adopt a simplified model for the re sonant Compton 
scattering (RCS. ILvutikov fe Gavriil 120061 ) of seed photons 
(assumed to have a pure blackbody spectral energy distri¬ 
bution) in NS magnetospheres, and s uccessfully appli ed to 
model spectra of several magnetars dRea et al-lboosl) . Al¬ 
though it is a ID, plane-parallel model, far from describing 
the non-trivial interaction between plasma and radiation in 
a realistic 3D magnetosphere, it effectively fits X-ray data. 
RCS models are parametrized by two parameters: the res¬ 
onant optical depth along a ray r res , which is related to 
the plasma density, and the therm al velocity of electrons in 
units of the speed of light, /3t (see ILvutikov fc Gavnilll200fi 
for more details). 

In Fig. [2] (left) we show the comparison between a black¬ 
body of kT = 1 keV, and the spectra predicted by the RCS 
model for /3t = 0.3 and different values of r res . We show the 
spectra as they would be seen when absorbed by the inter¬ 


stellar medium with Nh = 3 x 10 22 cm -2 , a typical value for 
magnetars in the galaxy. Since the interstellar absorption is 
important below ~ 1 keV, an efficient RCS could enhance 
the NSs detectability, because it promotes low energy pho¬ 
tonstojiigher energ ies. 

iRea et ahl f2008j) fit the X-ray spectra in quiescence of 
several magnetars. We found a correlation between B and 
the RCS parameters fir and r TBa , which can be reproduced 
by the following analytical fits: 


Pt = | 

0.001 

0.3 


f 0.001 

'Pres = 

{ B 


O 

T 

O 


if logR < 13 
if log B > 13 

if log B < 13 
if log B > 13 


(3) 

(4) 


We implement this analytical fit to modify the spectra 
of the randomly generated NSs in the population synthesis 
code. RCS does not change the original blackbody spectrum 
in NSs with low magnetic field, but it generates a hard tail 
in magnetars. In order to show the effect of the RCS in a 
typical PS realisation, in the middle panel of Fig. ?? we 
show the ratio of the RCS flux to the pure blackbody flux 
in the [0.1, 10] keV band, as a function of the blackbody 
flux. This figure clearly shows the main effect of the RCS 
correction: a flux enhancement by a factor of a few (up to 
a factor of 20) is applied to a significant number of sources, 
which has a visible imprint in the statistics of the number of 
detectable NSs. This effect is more evident for the dimmest 
sources: since RCS affects the high-energy tail, we can only 
obtain a very large enhancement for sources with a large Nh, 
and these are always dim. The right panel of Fig. [2] shows 
the same ratio as a function of the initial field strength. It 
shows that only for the highest fields, in the magnetar range, 
a significant flux enhancement in the [0.1-10] keV band is 
produced. 


3 RESULTS 

We begin by analyzing the solutions that best fit the radio¬ 
pulsar distribution (Tabic [I}. In Fig. [3] we compare the ac¬ 
cumulated number of detected X-ray sources N(S) with flux 
S > Sx s (hereafter log A-log S diagram) for each magne- 
tothermal evolution model. The sample of sources is not 
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Log S (erg cm 2 s ') 



Log Period (s) 

Figure 3. Log N - Log S diagrams (top) and period distributions 
( bottom ) for the thermally emitting X-ray pulsars for models A 
to F and the parameters listed in Table |T] In the bottom panel 
we show the histogram of sources with Sxmin > 10“ 14 erg s —1 
cm -2 . In colour lines we show the population synthesis results 
while the black line shows the observed sources. 

complete, since a full all-sky survey is not available for the 
present instruments, and many sources, in particular mag- 
netars, have been discovered while they are in outburst. We 
can assume that we are actually observing all of the few 
brightest sources, and we gradually miss more and more 
sources as the threshold flux decreases. From a quick in¬ 
spection of the log iV-log S diagrams, we can immediately 
conclude that none of the models of Table 1 predicts a num¬ 
ber of thermally emitting X-ray sources compatible with the 
observations. However, as shown in Fig. [T| there is a quite 
large region in the parameter space with similar statistical 
significance for the 25-value. We can then explore whether 
other parameter combinations (always within the region al¬ 
lowed by the radio-pulsar population analysis) can improve 
the log IV-log S results. Since we obtain a systematic lack of 
bright sources for all models, we have to modify the initial 
distribution to include more NSs born with high magnetic 
fields, which can be done if we move up-right in the dark 
blue region of Fig. [Tf 

We have proceeded as follows: first, we search for a pair 
of values of [j,b 0 , <jb 0 in the dark blue region of Fig. [l]than 
can reproduce the observed distribution of X-ray pulsars for 



Log S (erg cm 2 s ') 



Log Period (s) 

Figure 4. Same as Fig. [3] for the parameters listed in Table [2] 


fluxes above ss 10 -12 erg s~ 4 cm -2 (the other 3 parameters 
do not affect the X-ray population). Once this is done, we 
fix the values of these two parameters and repeat the min¬ 
imization procedure to find new values of fj,p 0 , ap 0 , and a 
that best fit the radio-pulsar distribution. In Table [2] and 
Fig. [4] (top) we show the results for the parameters that 
best reproduce the log N-\og S diagram for all models ex¬ 
cept model A, for which field decay was not included, and 
therefore the luminosity is not affected by the magnetic field. 
In other words, the standard cooling of a neutron star with¬ 
out including any additional heating by magnetic field de¬ 
cay cannot be reconciled with the X-ray data. For the other 
models, we find acceptable solutions with similar statistical 
significance for the radio-pulsar fits. The 23-value obtained 
with the 2D KS test shown in the last column is an average 
of 10 realizations, and we can see that they are compatible 
with the results of Table 1 within 1- or 2-a. 

Below 10 -12 erg s _1 cm~ 2 , the observed number of 
sources saturates due to the sensitivity limit of the present 
instruments, combined with the interstellar medium absorp¬ 
tion, which only allows us to observe dim sources if they 
are nearby. Thus, we conclude that the degeneracy shown 
in Fig. [I] where a large region of the parameter space was 
allowed by the radio-pulsar data, is severely reduced when 
the thermally emitting X-ray pulsars are included in the 
analysis. Interestingly, for different models we obtain sim- 
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Model 

^Bq 

log B [G] 

a Bo 

log B [G] 

[b] 

VPo 

[b] 

a 

n br 

[century -1 ] 

V 

V 

B 

13.34 

0.76 

0.30 

0.29 

0.44 

4.98 ±0.11 

3.1 

0.089 ± 0.009 

C 

13.36 

0.80 

0.29 

0.37 

0.44 

4.37 ±0.12 

4.0 

0.093 ±0.011 

D 

13.25 

0.77 

0.28 

0.20 

0.43 

3.07 ±0.11 

4.0 

0.089 ±0.010 

E 

13.23 

0.72 

0.32 

0.19 

0.44 

3.58 ±0.10 

2.5 

0.083 ± 0.008 

F 

13.15 

0.70 

0.29 

0.13 

0.44 

3.30 ±0.08 

5.0 

0.087 ±0.012 


Table 2. Optimal parameters for each magneto-thermal evolution model (defined in Table [TJ that fit the observed log AT - log S' distri¬ 
bution. The parameter 77 (see eq. 0 is also shown. 


ilar solutions for the initial magnetic field distribution, in 
the narrow region fis 0 = 13.15 — 13.35, ob 0 ~ 0.7 — 0.8 and 
a = 0.43 - 0.44. 

However, to reproduce the log IV-log S diagram is a nec¬ 
essary, but not sufficient condition to test a model. We now 
turn to explore in more detail the particular form of the 
period distributions. We remind again that the comparison 
is made with an incomplete observational data set, because 
the number of observed sources is only a fraction of the po¬ 
tentially detectable number of NSs. In Fig. [4] (bottom) we 
show spin period histograms comparing the observed sample 
(black solid line) with synthetic models. 

An important result can already be seen in these dis¬ 
tributions. On one hand, one needs a significant fraction of 
magnetars/high-B field NSs in the initial distribution to ac¬ 
count for the number of X-ray sources already present in 
our observed sample, but on the other hand, if that frac¬ 
tion is too large, the model over-predicts the number of 
X-ray pulsars with long periods, even assuming the exis- 
tence of a highly resistive layer in the inner crust of NSs 
(IPons. Vigano fc Real 120131 '). We can adjust the maximum 
period reached by a NS by varying the impurity parameter 
in the innermost part of the crust. Increasing the value of 
Qimp results in the decrease of the maximum period, but 
still not enough to successfully reproduce the observed up¬ 
per limit of 12 s and, at the same time, explain the observed 
number of very bright sources. The strong constraint result¬ 
ing from combining the need of a few very bright sources 
with the lack of isolated X-ray pulsars with periods longer 
than 12 s turns out to be a very restrictive condition. The 
next relevant question is whether or not unknown observa¬ 
tional selection effects are causing this discrepancy. 

To explore how the period distribution varies with flux, 
let us focus on model D (the other models are qualitatively 
similar). In Fig. [S] we compare the observed distribution to 
the synthetic data with a different cut-off in the absorbed 
flux (10 —14 , 10 -13 , and 10 -12 erg s -1 cm -2 ). An interesting 
trend is the apparent bimodal distribution of periods, similar 
to the observed distribution, although with a different scale 
and position. There is a partial lack of NSs with Pm l — 
2 s, as seen in the observations, despite the initial period 
distribution being a single Gaussian. However, even with the 
highest flux cut-off (10 -12 erg s -1 cm -2 ) there are always 
too many visible X-ray pulsars with long periods. Thus, this 
discrepancy cannot be attributed only to selection effects. 

Since it is impossible to a priori take into account (theo¬ 
retically) unknown selection effects in A'-rays, we propose in 
the following a phenomenological approach. A close inspec- 



Log Period (s) 


Figure 5. Same as Fig. [f] (bottom) for model D considering dif¬ 
ferent flux thresholds: 10 -12 (blue), 10 -13 (yellow) and 10 -14 
(red). All fluxes are in units of erg s —1 cm" 2 

tion of Fig. [4] shows that the observational sample appears 
to be complete (given all the assumptions made in paper and 
within statistical fluctuations) for fluxes Sx 13 > 3 x 10 -12 
erg s -1 cm -2 . We have also empirically found that the ratio 
of the number of observed sources to the number of poten¬ 
tially visible sources in the synthetic sample scales linearly 
with the flux in the range 10 -13 — 10 -12 erg s -1 cm -2 . If we 
assume that the probability to detect an X-ray pulsar with 
a given flux is simply 

Pobs = min {pS-n, 1} (5) 

where S_n is the absorbed flux in units of 10 -11 erg s -1 
cm -2 , we can obtain for each model the value of p that 
reproduces the log N — log Sx 13 diagram. It turns out that 
rj ~ 3 — 5 for all models (eight column in Tabic [2]). This 
naive approach is useful to estimate the number of missed 
sources by uncontrolled selection effects. 

We can now generate a random observed sample from 
the synthetically generated sample by imposing that the 
probability to detect the thermal emission of a pulsar with a 
given flux is given by the above expression. The correspond¬ 
ing log JV-log S diagrams and period distributions after such 
selection is applied are shown in Fig. [6] Although most pul¬ 
sars with high periods are now not observed , we still find 
a non-negligible number of bright enough sources with long 
periods. Our predictions range from 4 to 20 X-ray pulsars to 
be observed with P > 12 s, depending on the model. Even 
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Log S (erg cm 2 s ') 



11 12 13 14 15 16 

Log B 0 (G) 



Log Period (s) 

Figure 6. Same as Fig. [T] but applying the phenomenological 
detection probability formula (eg. l5l. 

considering that we are dealing with low-number statistics, 
this appears to be a real problem. 

3.1 Beyond Gaussian distributions of initial 
magnetic fields. 

All previous results (and most studies in the literature) rely 
on the assumption that the initial magnetic field distribution 
is log-normal. As the initial magnetic field increases, old NSs 
reach longer periods and therefore the narrower beam size 
produces a strong selection of radio-pulsars. Therefore, the 
low field part of the distribution is what we can practically 
constrain using the radio-pulsar population, while the high 
field part of the distribution remains unconstrained, because 
the distribution of high field radio-pulsars is highly degener¬ 
ate with other parameters (beaming angle, radio-luminosity 
correlation, etc.). We have seen that the optimal parameters 
for radio-pulsars underpredicts the number of X-ray sources 
with high fluxes. We hence need to increase the width of 
the Gaussian to generate more bright NSs with high mag¬ 
netic fields, but the problem is that the Gaussian becomes 
too wide, and we run into an overprediction of X-ray pul¬ 
sars with periods P > 12s, which has not been observed. 
The easiest, and perhaps more obvious, way out is to re¬ 
lax the assumption of an initial log-normal magnetic field 


Figure 7. Bo distributions for model D that best fit the radio¬ 
pulsar (green/solid line) and X-ray pulsar (brown/dot-dashed) 
populations. In red/three-dot dashed and blue/dashed we plot 
the results for truncated and bimodal distributions discussed in 
Section |3.II respectively. 

distribution. For example, if we consider model D, we can 
visualize the Bo distribution functions for the best fit to the 
radio-pulsar population (green/solid) and X-ray populations 
(brown/dot-dashed) in Fig. [7] The excess of long period X- 
ray pulsars is originated by the high Bo tail of the latter. 
We now discuss two possibilities, to reconcile the synthetic 
models with the data. 

3.1.1 A truncated magnetic field distribution. 

If the origin of the neutron star magnetic field is attributed 
to magneto-hydrodynamic instabilities (e.g., any dynamo 
process related to rotation or convection), a nonlinear sat¬ 
uration is expected when the amplified magnetic field be¬ 
comes dynamically relevant to suppress the instability. This 
may happen at typical fields of t he order of 10 15 G, d e pend¬ 
ing on the particular mecha nis m llObergaulinger et alJl2009l : 
lObergaulinger. Janka fc Alov! l2014h . Thus. we can simply 

introduce a cut-off value Bo max to the log-normal distribu¬ 
tion (red/three-dot dashed curve in Fig. [TJ • In Fig. [8] we can 
see that the log A-log S diagram for the X-ray population 
can also be successfully reproduced by a truncated distribu¬ 
tion with Bomax = 5 x 10 14 G and increased p,B 0 and <tb 0 
(Table O. The period distribution is significantly modified, 
and now there are no visible sources with P > 20 s, as shown 
in Fig. [8] 

3.1.2 A second population of magnetised NSs. 

The alternative possibility that we have explored is the ex¬ 
istence of a second population of magnetised NSs, whose 
origin is different from that of standard radio-pulsars. The 
reason for this second population could be attributed to two 
different evolutionary paths of massive stars if they are iso¬ 
lated or in binary systems. Since about 80% of massive stars 
are in binaries, and assuming that half of them go through 
a common envelope phase or some evolutionary stage that 
produces a magnetized stellar core, one could expect about 
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Model 

^ Bq 

log B [G] 

a Bo 

log B [G] 

[b] 

°P 0 

W 

a 

n br 

[century -1 ] 

V 

V 

D 

13.75 

1.15 

0.44 

0.26 

0.44 

6.81 ±0.20 

3.1 

0.120 ±0.012 

E 

13.58 

0.90 

0.33 

0.21 

0.44 

5.31 ±0.19 

5.0 

0.091 ±0.009 

F 

13.33 

0.83 

0.30 

0.19 

0.43 

3.24 ± 0.09 

5.0 

0.110 ±0.008 


Table 3. Optimal parameters for models D, E and F with a truncated distribution of initial magnetic fields (Bonwi 


5 x 10 14 G). 




Figure 8. Same as Fig.[6]assuming a log-normal distribution with 
cutoff at Bo m ax = 5 X 10 14 G, and the parameters of Table m 


Model 

NS1 

(%) 

NS2 

(%) 

7T.br 

[century -1 ] 

V 

V 

D 

60 

40 

4.34 ±0.13 

6.3 

0.087 ±0.007 

E 

65 

35 

4.80 ±0.10 

2.0 

0.082 ± 0.006 

F 

70 

30 

3.78 ±0.08 

3.2 

0.081 ± 0.005 


Table 4. Results for models D, E and F with a bimodal dis¬ 
tribution. We show the relative weight of the first population 
(NS 1, same parameters as in Table [T} and second population 
of highly-magnetized stars (NS2, flat distribution in the range 
log S 0 G [13.5,14.5] ). 


40% of the total pulsar population being born with a bi¬ 
nary massive progenitors, poss ibly leading to the forma- 
tion of more magnetized NSs (ISpruitl 120081 ; lLangerl 120121 : 
IClark et, al.ll2014h . We have investigated what happens if, 
in addition to the population generated with the param¬ 
eters of Table [3] by fitting only the radio pulsar popula¬ 
tion (hereafter denoted by NS1), we add a second popu¬ 
lation (denoted by NS2), uniformly distributed in the range 
log Bo G [13.5, 14.5] and normalized to account for the total 
number of NSs observed (blue/dashed curve in Fig. Q. The 
results are shown in Fig. [9] where again we can observe that 
both the log IV-log S diagram and the period distribution 
are satisfactorily reproduced. We have verified that this sec¬ 
ond component does not modify the radio pulsar population, 
since their magnetic fields are so large that the vast major¬ 
ity of NS2 sources remains undetectable as radio-pulsars. 
The total birthrate, however, will increase accordingly, still 
marginally compatible with the observations. We also note 
that one can obtain similar results by assuming other dis¬ 
tributions for this second population, for example a narrow 
log-normal distribution centered in « 1 — 3 x 10 14 G. 


3.2 The upper limit to the birth rate of very 
strongly magnetized NSs. 

Before concluding our analysis on the period distribution, a 
simple calculation can stress once more the conflict between 
having a significant fraction of NSs with magnetic fields at 
birth above a certain value, and the lack of observed isolated 
NSs with periods > 12 s. To quantify this incompatibility, we 
begin by obtaining the probability of observing a star, born 
with a given initial magnetic field Bo, as an X-ray pulsar 
with a period longer than 12 s. The results for three different 
values of the initial magnetic field are given in Table [S] for 
the representative model D (similar values are obtained for 
other models). In this calculation, we assumed a uniform 
distribution of ages in the range 0-560 Myrs (with constant 
birth rate), and we generate a large number of NSs with the 
same magnetic field. Then, we can obtain the probability as 
the fraction of synthetic visible sources with P > 12 s to the 
total number of stars generated. 

Next we address the following question: what is the 
maximum number of stars born with Bo compatible (within 
statistical fluctuations) with the absence of observed neu¬ 
tron stars with periods longer than 12 s? This number can 
be estimated taking into account that if Nb 0 is the number 
of stars born with a given magnetic field, the probability 
of observing none of them with periods longer than 12 s is 
given by e _JVs o p , according to the Poisson distribution, and 
the probability that we will observe at least one star with 
P > 12 s is 1 - e~ NB o p . 
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Log Period (s) 


Figure 9. Same as Fig. [6] with a bimodal distribution of Bo 
described in §3.1.2. 


We can use this value to reject the null hypothesis in the 
following sense: for a given model, we can calculate the value 
of Nb 0 that corresponds to a probability of 1 — e~ NB o p = 
0 .01, and we know that in 99% of the realizations with Nb 0 
generated stars, we will not detect any source with P > 
12 s. Therefore, we can conclude that if a model predicts 
^ Nb 0 stars, the fact that none has been observed rejects the 
model at the 99% confidence level (statistical fluctuations 
may result in a no-detection only in ^ 1% of the cases). 

In Fig. [10] we plot the probability of no-detection of NSs 
with periods longer than 12 s against the birthrate, for dif¬ 
ferent values of Bo . This can be interpreted as an upper limit 
to the birthrate of stars with initial magnetic field > Bo, for 
a given confidence level. For a confidence level of 99% and a 
magnetic field of 10 15 G we obtain Nb 0 = 112000. This value 
corresponds to a birthrate of ~ 0.02 objects/century which 
means, by comparison with the total galactic NS birthrate 
~ 3 NSs/century, that less than 0.7% can be born with mag¬ 
netic fields > 10 15 G. Of course, this fraction is even smaller 
if we reduce our confidence level. The maximum birthrate 
for stars with initial magnetic field greater than 5 x 10 14 G 
can be 3 times larger than for Bo = 10 15 G for the same 
confidence level. This is satisfied by both alternative mod¬ 
els described in the two subsections above, which cannot be 
discriminated with current data. 


Bo [G] 

P 

5 x 10 14 

1.3 x lO” 5 

8 x 10 14 

3.1 x 1CT 5 

1 x 10 15 

4.1 x 10-® 


Table 5. Probability of detection of stars whose period is longer 
than 12 s for samples with fixed initial magnetic field Bq. 



0.00 0.01 0.02 0.03 0.04 0.05 

BR Mognetor (NS/century) 


Figure 10. Probability of no-detection of NSs with P > 12 s 
as a function of the birthrate, for three different initial magnetic 
fields: So = 1 X 10 15 G (blue squares), 8x 10 14 G (green triangles) 
and 5 X 10 14 G (red diamonds). The symbols show results of 
simulations with 20 million NSs and the dashed lines correspond 
to the theoretical Poissonian value e~ B o p . 


3.3 P distribution. 

To complete our discussion, we comment on our results for 
the distribution of P. In Fig.[lT]we show the histograms for 
the models considered, which roughly agree with the obser¬ 
vations. However, all models and initial field distributions 
(log-normal, truncated log-normal, bimodal) show a similar 
trend: synthetic samples display an excess of sources with 


with P>10 11 ss 1 . This is due to the fact that even 
the strongest magnetic fields [Bo > 10 15 G) hardly achieve 
P ~ 10 -10 s s _1 . 

The observed sources with very high P are mag- 
netars and it is likely that modelling their spin evo- 
lutio n with th e sim p le magnetospheric torque formu¬ 
l ae (ISpi tkovskyl 2006j; Beskin ^ Istomin fe Philippovl 120131 : 
IPhilippov, Tchekhovskov fe Lil 120131) is an oversimplifica¬ 
tion. For instance, strong particle winds can_b e non- 

negli gible in the early ~ 1000 yr of their lives (jTong et al.l 
120131) when these sources undergo frequent bursts and flares. 
As a result of additional torques, the P can increase by an 
order of magnitude for the same magnetic field. This would 
imply that the magnetic field estimate by the classical mag- 
netodipolar formula is a factor of a few higher than the real 
value. This effect is only expected to work during a short 
time and does not change the results of the analysis of the 
whole NS population, but it may contribute to the anoma¬ 
lous high values of P of the few youngest magnetars. 
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Log P 




Log P 

Figure 11. Observed and predicted P distributions. Top: mod¬ 
els with log-normal distributions. Middle: models with truncated 
distributions discussed in section [3.1.II Bottom: models with bi- 
modal distributions discussed in Section 13.1.21 In colour lines 
we show the PS results while the black line shows the observed 
sources. 

4 SUMMARY 

We have extended previous works on population synthesis 
of isolated neutron stars to study the combination of initial 
conditions and microphysical parameters required to explain 
the radio-pulsars and thermally emitting A'-ray pulsar pop¬ 
ulations within the same evolutionary model. 


Our first finding is that the log A-log S diagram of X- 
ray pulsars can be well fitted by several models that, at the 
same time, reproduce the radio-pulsar distribution (because 
the parameter space of initial magnetic field distribution of 
radio-pulsars has a large region with the same statistical sig¬ 
nificance) . Models with iron envelopes require a higher mean 
value of the initial magnetic field and a wider distribution 
than those with light elements. 

Our second relevant result is that the obtained sample 
of observed sources with Stf s > 3 x 1CP 12 erg s -1 cnY 2 is 
nearly complete, and we estimate that, with more sensitive 
instruments, the thermal emission of about one hundred of 
NSs with fluxes Stf s > 10 -13 erg s _1 cm -2 is potentially 
detectable. This number can be increased up to 500 sources 
if the sensitivity is good enough to cover the whole sky with 
Sf s > 10~ 14 erg s’ 1 cm" 2 . 

However, when examining in detail the period distribu¬ 
tion all models with log-normal distributions of the initial 
magnetic field fail to explain why there are no X-ray pul¬ 
sars with periods longer than 12 seconds. This sharp upper 
limit cannot be reproduced, even with models with a very 
high resistivity in the crust/core interface, if the initial field 
distribution is assumed to be log-normal. The problem is 
solved if other distributions are invoked, such as a truncated 
log-normal distribution, or a binormal distribution with two 
distinct populations. The most interesting feature, common 
to the two alternative cases considered, is that both cases 
require that magnetars cannot be born with a very high 
magnetic field. For a confidence level of 99%, we obtain an 
upper limit to the birth rate of NSs with Bo > 10 15 G of 
< 0.02 objects/century, which represents a tiny fraction of 
the population (less than 1%). 

The existence of an upper limit to the magnetic field 
at birth in the range B o = 5-8x 10 14 G, combined with 
fast magnetic field dissipation and observational selection 
effects, explains the lack of isolated pulsars with long peri¬ 
ods. Although this has a drawback, the lack of a few pul¬ 
sars with very high values of P > 10 —11 , note that we have 
assumed that NSs spindown only by the rotating magneto- 
spheric torque. In very young, active magnetars, losses of 
angular momentum carried by p articles (wind ) can compete 
with the magnetospheric torque llTong et alil2013l l , increas¬ 
ing the value of P for a given magnetic field. A more detailed 
study of these phenomena, and other modifications of the 
spin-down evolution is reserved to future works. 
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